//Features as per Andres' WACV paper

inline
opt_feat::opt_feat(const std::string in_path,
		   const std::string in_actionNames,  
		   const uword in_col, 
		   const uword in_row
)
:path(in_path), actionNames(in_actionNames), col(in_col), row(in_row)
{
  actions.load( actionNames );  //actions.print("All actions");
  ismultiAction = false; // per default is false
  dim = 14;
}



inline
void
opt_feat::features_per_action_training( field<string>peo_train ) 
{
  std::vector < vec > feat_all_videos_action_i;
  int n_actions = actions.n_rows;
  all_actions_matrix.resize(n_actions);
  //cout << "size " << all_actions_matrix.size() << endl;
  
  /*for (uword b=0;  b< all_actions_matrix.size(); ++b)
   *    {
   *      cout << all_actions_matrix.at(b).size() << " ";
}
*/cout << endl;


int n_peo_tr =  peo_train.n_rows;

for (int pe = 0; pe< n_peo_tr; ++pe)
{
  
  all_actions_matrix.clear();
  all_actions_matrix.resize(n_actions);
  
  //esto cambia
  std::stringstream ss_video_name;
  ss_video_name << path <<  peo_train (pe) << "_Brownie.avi";
  
  std::stringstream ss_label_name;
  ss_label_name << path <<  peo_train (pe) << "_unique_verbs_labels.dat";
  
  std::stringstream ss_range_frames;
  ss_range_frames << path <<  peo_train (pe) << "_start_end.dat";
  
  
  cout << ss_video_name.str() << endl;
  
  feature_video( ss_video_name.str(), ss_label_name.str(), ss_range_frames.str() ) ; //all_actions_matrix is calculated inside this method
  
  
  
  // 	cout << "size " << all_actions_matrix.size() << endl;
  // 	for (uword b=0;  b< all_actions_matrix.size(); ++b)
  // 	{
  // 	  cout << all_actions_matrix.at(b).size() << " ";
  // 	  
  // 	}
  // 	
  // 	cout << endl;
  
  for (int act=0; act<n_actions; ++act)
  {
    feat_all_videos_action_i = all_actions_matrix.at(act);
    //cout << feat_all_videos_action_i.size() << endl;
    //cout << feat_all_videos_action_i.at(0).size() << endl;
    
    
    mat mat_features;
    if (feat_all_videos_action_i.size()>0)
    {
      mat_features.zeros(dim,feat_all_videos_action_i.size());
      //cout << mat_features.n_rows << " " << mat_features.n_cols << endl;
      
      for (uword i = 0; i < feat_all_videos_action_i.size(); ++i)
      {
	mat_features.col(i) = feat_all_videos_action_i.at(i);
	//cout << i << " " ;
	
      }
    }
    else
    {
      mat_features.zeros(dim,0);
    }
    std::stringstream tmp_ss4;
    tmp_ss4 << "./features/train/feature_vectors_dim" << dim << "_" << peo_train (pe) << "_" << actions(act);
    
    mat_features.save( tmp_ss4.str() );
    
  }
  
}

}


///Create universal GMM Model 

inline
void
opt_feat::create_universal_gmm(field<string>peo_train, int in_Ncent, int run)
{ 
  N_cent = in_Ncent;
  cout << "Calculating Universal GMM " << endl;
  cout << "# clusters: " << N_cent << endl;
  
  mat uni_features;
  
  
  for (uword act = 0 ; act < actions.n_rows;  ++act) {
    
    mat mat_features_tmp;
    mat mat_features;
    
    for (uword pe=0; pe<peo_train.n_rows; ++pe)    {
      mat mat_features_person_i;
      
      
      std::stringstream tmp_ss4;
      tmp_ss4 << "./features/train/feature_vectors_dim" << dim << "_" << peo_train (pe) << "_" << actions(act);
      //tmp_ss4 << "./features/train/feature_vectors_" << peo_train (pe) << "_" << actions(act);
      //cout << "Doing for " << tmp_ss4.str() << endl;
      mat_features_person_i.load( tmp_ss4.str() );
      
      if ( mat_features_person_i.n_cols>0 )
      {
      mat_features_tmp	 = join_rows( mat_features_tmp, mat_features_person_i );
      }
      else 
      {
	cout << "# vectors = 0 in " << tmp_ss4.str() << endl;
      }
      
     //mat_features_tmp	 = join_rows( mat_features_tmp, mat_features_person_i );
     //cout << "rows & cols " << mat_features.n_rows << " & " << mat_features.n_cols << endl;
      
    }
    
    const uword N_max = 100000;  // maximum number of vectors per action to create universal GMM
    
    if (mat_features_tmp.n_cols > N_max)
    {
      ivec tmp1 = randi( N_max, distr_param(0,mat_features_tmp.n_cols-1) );
      ivec tmp2 = unique(tmp1);
      uvec my_indices = conv_to<uvec>::from(tmp2);
      mat_features = mat_features_tmp.cols(my_indices);  // extract a subset of the columns
      
    }
    else
    {
      mat_features = mat_features_tmp;
      
    }
    
    
      if ( mat_features.n_cols>0 )
      {
	uni_features	 = join_rows( uni_features, mat_features );
	
      }
      else 
      {
	cout << "# vectors = 0 in uni_features" << endl;
	
      }
    
    //uni_features =  join_rows( uni_features, mat_features );
    mat_features_tmp.reset();
    mat_features.reset();
    
    
  }
  
  
  cout << "r&c "<<  uni_features.n_rows << " & " << uni_features.n_cols << endl;
  bool is_finite = uni_features.is_finite();
 
    if (!is_finite )
    {
      cout << "is_finite?? " << is_finite << endl;
      cout << uni_features.n_rows << " " << uni_features.n_cols << endl;
      getchar();
    
    }
  
  
  jo_gmm_diag gmm_model;
  jo_gmm_diag bg_model;
  
  
  
  bool status_em = false;
  int rep_em=0;
  
  while (!status_em)
  {
    bool status_kmeans = false;
    //int rep_km = 0;
    
    while (!status_kmeans)
    {
      arma_rng::set_seed_random();
      
      status_kmeans = gmm_model.learn(uni_features, N_cent, jo_eucl_dist, jo_random_subset, 10, 0, 1e-5, false, false);   //Only Kmeans
      bg_model = gmm_model;
      //rep_km++;
    }
    
    
    status_em = gmm_model.learn(uni_features, N_cent, jo_eucl_dist, jo_keep_existing, 0, 10, 1e-5, true, false);   
    rep_em++;
    
    if (rep_em==9)
    {
      status_em = true;
      gmm_model = bg_model;
      
    }
    
  }
  
  
  cout <<"EM was repeated " << rep_em << endl;
  
  std::stringstream tmp_ss5;
  tmp_ss5 << "./run"<< run <<"/gmm_models/UniversalGMM_Ng" << N_cent << "_dim" <<dim; 
  cout << "Saving GMM in " << tmp_ss5.str() << endl;
  gmm_model.save( tmp_ss5.str() );
  cout << endl;
  
  
}

///Create GMM Model per Action  

inline
void
opt_feat::create_gmm_action(field<string>peo_train, int in_Ncent, int run)
{ 
  N_cent = in_Ncent;
  cout << "# clusters: " << N_cent << endl;
   
  std::stringstream tmp_ss5;
  tmp_ss5 << "./run"<< run <<"/gmm_models/UniversalGMM_Ng" << N_cent << "_dim" <<dim; 
  jo_gmm_diag universal_gmm_model;
  universal_gmm_model.load( tmp_ss5.str() );
  
  
  for (uword act = 0 ; act < actions.n_rows;  ++act) {
    
    cout << "Calculating GMM for action " << actions(act) << endl;
    mat mat_features_tmp;
    mat mat_features;
    
    for (uword pe=0; pe<peo_train.n_rows; ++pe)    {
      mat mat_features_person_i;
      
      
      std::stringstream tmp_ss4;
      tmp_ss4 << "./features/train/feature_vectors_dim" << dim << "_" << peo_train (pe) << "_" << actions(act);
      //tmp_ss4 << "./features/train/feature_vectors_" << peo_train (pe) << "_" << actions(act);
      //cout << "Doing for " << tmp_ss4.str() << endl;
      mat_features_person_i.load( tmp_ss4.str() );
      
      if ( mat_features_person_i.n_cols>0 )
      {
      mat_features_tmp	 = join_rows( mat_features_tmp, mat_features_person_i );
      }
      else 
      {
	cout << "# vectors = 0 in " << tmp_ss4.str() << endl;
      }
      
    }
    
    const uword N_max = 200000;  // maximum number of vectors
    if (mat_features_tmp.n_cols > N_max)
    {
      ivec tmp1 = randi( N_max, distr_param(0,mat_features_tmp.n_cols-1) );
      ivec tmp2 = unique(tmp1);
      uvec my_indices = conv_to<uvec>::from(tmp2);
      mat_features = mat_features_tmp.cols(my_indices);  // extract a subset of the columns
    }
    else
    {
      mat_features = mat_features_tmp;
      
    }
    
    mat_features_tmp.reset();
    
    
    cout << "r&c "<<  mat_features.n_rows << " & " << mat_features.n_cols << endl;
    bool is_finite = mat_features.is_finite();
    
    if (!is_finite)
    {
      cout << "is_finite?? " << is_finite << endl;
      cout << mat_features.n_rows << " " << mat_features.n_cols << endl;
   
      getchar();
    
    }
    
    
    jo_gmm_diag gmm_model = universal_gmm_model;
    
    
    
    bool status_em = false;
    //int rep_em=0;
    
    int km_iter = 0;
    int em_iter = 10;
    int var_floor = 1e-5;
    bool print_mode = true;
    bool means_only = true;
    
    
    
    status_em = gmm_model.learn(mat_features, N_cent, jo_eucl_dist, jo_keep_existing, km_iter, em_iter, var_floor, print_mode, means_only);   //only_means
    
    if (status_em==false)
    {
      cout << "Saving Matrix. See what's the problem " << endl;
      mat_features.save("problem_here_2.dat", raw_ascii );
      cout << "EM failed" << endl;
      
    }
      
      //rep_em++;

    
    
    //cout <<"EM was repeated " << rep_em << endl;
        
    std::stringstream tmp_ss5;
    tmp_ss5 << "./run"<< run <<"/gmm_models/Ng" << N_cent << "_" << actions(act)<< "_dim" <<dim; 
    cout << "Saving GMM in " << tmp_ss5.str() << endl;
    gmm_model.save( tmp_ss5.str() );
    cout << endl;
    
  }
}


// //*******************************************************************************************
// //***************************Feature per video: Training and Testing ************************

inline 
void
opt_feat::feature_video(std::string one_video, std::string sLabel, std::string sRange_frames)
{
  
  std::vector < vec > feat_frame;
  
  vec labels; 
  labels.load( sLabel );
  //cout << " " << labels.n_elem << endl;
  
  
  //labels.t().print("labels");
  
  rowvec range_frame;
  range_frame.load( sRange_frames );
  //range_frame.print("ini end");
  int ini_frame = range_frame(0);
  //int end_frame = range_frame(1) - range_frame(0);
  int con_fr_tr=0;
  int con_fr_te=0;
  
  
  
  
  cv::VideoCapture capVideo(one_video);
  //cout << one_video << endl;
  //double fps = capVideo.get(CV_CAP_PROP_FPS); //get the frames per seconds of the video
  //cout << "Frame per seconds : " << fps << endl;
  
  //cv::namedWindow("MyVideo",CV_WINDOW_AUTOSIZE); //create a window called "MyVideo"
  
  //double frmcount = capVideo.get(CV_CAP_PROP_FRAME_COUNT);
  capVideo.set(CV_CAP_PROP_POS_FRAMES,ini_frame);
  
  //cout << "# of frames is: " << frmcount << endl;
  
  
  
  if( !capVideo.isOpened() )
  {
    cout << "Video couldn't be opened" << endl;
    return;
  }
  
  cv::Mat prevgray, gray, flow, cflow, frame, prevflow;
  cv::Mat ixMat, iyMat, ixxMat, iyyMat;
  cv::Mat flow_xy[2], mag, ang;
  float Gten, Sten;
  //cv::namedWindow("My Video", 1);
  //running_stat_vec<vec> stats_video(true);
  //cout << "Frame: ";
  
  //Saving each four frames
  int n_frames = labels.n_elem;
  //cout << "n_frames " << n_frames << endl;
  //int len = ceil(n_frames/4); 
  int len = ceil( (n_frames - 1)/2 );
  //int len = n_frames - 2; //0 and 1 not used
  //len=len-1; //0 and 1 not used
  //int len = n_frames - 2 ;//0 & 1 are not used
  //cout << "Len " << len << endl;
  test_features_arma.set_size( len );
  lab_feature_vectors.set_size( len );
  
  
  
  //cout << "end_frame " << end_frame << endl;
  int wfr;
  for(int fr=0; fr<n_frames; fr++){
    
    //cout << fr << " " ;
    wfr = fr;
    bool bSuccess = capVideo.read(frame); // read a new frame from video
    
    if (!bSuccess) //if not success, break loop
    {
      cout << "Cannot read the frame from video file" << endl;
      break;
    }
    
    
    cv::resize(frame,frame, cv::Size(col,row) );
    cv::cvtColor(frame, gray, CV_BGR2GRAY);
    
    
    
    
    if( prevgray.data )
    {
      
      cv::calcOpticalFlowFarneback(prevgray, 
				   gray, 
				   flow, 
				   0.5, //pyr_scale
				   3,   //levels
				   9,   //winsize
				   1,   //iterations
				   5,   //poly_n
				   1.1, //poly_sigma
				   0);  //flags
      
      
      //optical flow
      cv::split(flow, flow_xy);
      cv::cartToPolar(flow_xy[0], flow_xy[1], mag, ang, true);
      
      
      cv::Sobel(gray, ixMat, CV_32F, 1, 0, 1);
      cv::Sobel(gray, iyMat, CV_32F, 0, 1, 1);
      cv::Sobel(gray, ixxMat, CV_32F, 2, 0, 1);
      cv::Sobel(gray, iyyMat, CV_32F, 0, 2, 1);
      
      float  ux = 0, uy = 0, vx = 0,  vy = 0;
      float u, v;
      
      if( prevflow.data )
      {
	
	for (uword x = 0 ; x < col ; ++x ){
	  for (uword y = 0 ; y < row ; ++y ) {
	    
	    vec features_one_pixel(dim);
	    u = flow.at<cv::Vec2f>(y, x)[0];
	    v = flow.at<cv::Vec2f>(y, x)[1];
	    
	    //cout << "x= " << x << " - y= " << y << endl;
	    // x grad
	    //cout << " x y grad" << endl;
	    float ix = ixMat.at<float>(y, x);
	    //cout << " y grad" << endl;
	    float iy = iyMat.at<float>(y, x);
	    
	    // grad direction &  grad magnitude
	    //cout << "grad direction &  grad magnitude" << endl;
	    float gd = std::atan2(std::abs(iy), std::abs(ix));
	    float gm = std::sqrt(ix * ix + iy * iy);
	    
	    // x second grad
	    //cout << "x y  second grad " << endl;
	    float ixx = ixxMat.at<float>(y, x);
	    // y second grad
	    float iyy = iyyMat.at<float>(y, x);
	    
	    //du/dt
	    float ut = u - prevflow.at<cv::Vec2f>(y, x)[0];
	    // dv/dt
	    float vt = v - prevflow.at<cv::Vec2f>(y, x)[1];
	    
	    //// divergence &  vorticity
	    //cout << "divergence &  vorticity" << endl;
	    if (x>0 && y>0 )
	    {
	      ux = u - flow.at<cv::Vec2f>(y, x - 1)[0];
	      uy = u - flow.at<cv::Vec2f>(y - 1, x)[0];
	      vx = v - flow.at<cv::Vec2f>(y, x - 1)[1];
	      vy = v - flow.at<cv::Vec2f>(y - 1, x)[1];
	    }
	    
	    
	    //14 features
	    features_one_pixel  << x << y << abs(ix) << abs(iy) << abs(ixx) 
	    << abs(iyy) << gm << gd <<  u << v << abs(ut) 
	    << abs(vt) << (ux + vy)  << (vx - uy);

	    
	    double is_zero = accu(abs(features_one_pixel));

	    if (is_finite( features_one_pixel ) && is_zero!=0 )
	    {
	    if (gm>40 )  // Empirically set to 40
	    {
	      
	      /*frame.at<cv::Vec3b>(y,x)[0] = 0;
	       *			       frame.at<cv::Vec3b>(y,x)[1] = 0;
	       *			       frame.at<cv::Vec3b>(y,x)[2] = 255;
	       */
	      feat_frame.push_back(features_one_pixel);
	    }
	    
	    }
	    else
	    {
	      //bool is_fin = is_finite( features_one_pixel );
	      //bool is_zer = is_zero!=0;
	      //cout << is_fin << " " << is_zer << endl; 
	      cout << "It's not FINITE or All elements are zero. Vector discarded " << endl;
	      cout << fr << endl;
	      features_one_pixel.t().print("Features Current Pixel: ");
	      
	    }
	    
	  }
	  
	}
	
	
	//cout << "Adding to Matrix" << endl;
	//cout << labels(fr) << endl;
	uword pos = labels(fr) - 1;
	if (ismultiAction)
	{
	  //cout << fr << " " ;
	  
	  
	  if ( (fr % 2 ) == 0 ) 
	  {
	    mat feat_frame_i( dim,feat_frame.size() );
	    for (uword l = 0; l<feat_frame.size(); ++l )
	    {
	      feat_frame_i.col(l) = feat_frame.at(l);
	      
	    }
	    test_features_arma(con_fr_te) = feat_frame_i;
	    lab_feature_vectors(con_fr_te) = pos;
	    
	    con_fr_te++;
	    //cout << con_fr_te << " ";
	  }
	  
	  
	}
	
	else{  
	  if ( (fr % 4 ) == 0 ) 
	  {
	    //feat_all_videos_action_i.insert(feat_all_videos_action_i.end(), feat_frame.begin(), feat_frame.end()  );
	    //cout << "before " << all_actions_matrix.at(pos).size() << " ";
	    all_actions_matrix.at(pos).insert(all_actions_matrix.at(pos).end(), feat_frame.begin(), feat_frame.end()  );
	    //cout << "after " <<all_actions_matrix.at(pos).size() << endl;
	    //getchar();
	    //all_actions_matrix.at(pos) = feat_frame;
	    
	  }
	  con_fr_tr++;
	}
	
	
      }
      
    }
    
    if(cv::waitKey(30)>=0)
      break;
    
    
    std::swap(prevgray, gray);
    std::swap(prevflow, flow);
    
    
    //Put frame number in each frame
    
    // 	  string text;
    //          stringstream tostring;
    //          tostring << fr << " " << actions( labels( fr ) - 1) ;
    //          text = tostring.str();
    //          int fontFace = cv::FONT_HERSHEY_PLAIN;
    //          double fontScale = 1.5;
    //          int thickness = 1;  
    //          cv::Point textOrg(20, 20);
    //          cv::putText(frame, text, textOrg, fontFace, fontScale, cv::Scalar::all(255), thickness,8);
    
    //cv::namedWindow("color", CV_WINDOW_NORMAL);
    //cv::setWindowProperty("color", CV_WND_PROP_FULLSCREEN, CV_WINDOW_FULLSCREEN);
    
    
    //cv::namedWindow("color",CV_WINDOW_FULLSCREEN);
    //cv::imshow("color", frame);    
    //cv::waitKey(1);
    
    feat_frame.clear();
  }
  
  //cout << " press a key " << endl;
  // cout << "con_fr_te " << con_fr_te << endl;
  // cout << "n_frames " << n_frames << endl;
  // cout << "len " << len << endl;
  // cout << "fr " << wfr << endl;
  
  //cout << "Descomentar. Press a key " << endl;
  //getchar();
  
}

// //******************TESTING MULTI ACTIONS***********************************
// //***************************************************************************
// 
// ///change this one :)
// 
// //load videos and labels list
// 
// 
inline 
void 
opt_feat::feature_multi_action( field<string> peo_test )
{
  ismultiAction = true; // per default is false
  
  int n_peo_test =  peo_test.n_rows;
  
  for (int pe = 0; pe <n_peo_test; ++pe ){ //videos.n_rows
    
    
    std::stringstream ss_video_name;
    ss_video_name << path <<  peo_test(pe) << "_Brownie.avi";
    
    std::stringstream ss_label_name;
    ss_label_name << path <<  peo_test(pe) << "_unique_verbs_labels.dat";
    
    std::stringstream ss_range_frames;
    ss_range_frames << path <<  peo_test (pe) << "_start_end.dat";
    
    
    cout << ss_video_name.str() << endl;
    feature_video(ss_video_name.str(), ss_label_name.str(), ss_range_frames.str() ) ;  //test_features_arma is calculated inside this method
    
    
    std::stringstream tmp_ss4;
    tmp_ss4 << "./features/multi_test/feat_dim" << dim << "_" << peo_test(pe);  
    
    std::stringstream tmp_vec_lab;
    tmp_vec_lab << "./features/multi_test/lab_"<< peo_test(pe);  
    
    
    //cout << tmp_ss4.str() << endl;
    //getchar();
    cout << "Saving " << endl;
    test_features_arma.save( tmp_ss4.str() );
    lab_feature_vectors.save( tmp_vec_lab.str(), raw_ascii );
    //cout << "Press a key" << endl;
    //getchar();
    
    
  }
  
}

inline 
void 
opt_feat::gmm_multi_action( field<string> peo_test, int in_Ncent, int L,  int run) 
{
  N_cent = in_Ncent;
  double ave_acc=0;
  
  field<mat> test_features_arma; // Features for frame i are in  arow of this field
  
  vec real_labels;
  
  cout << "Testing for GMM with " << N_cent << " centroids" << endl;
  
  
  vec likelihood_actions;
  double universal_likelihood;
  vec likelihood_label;
  
  int n_peo_test =  peo_test.n_rows;
  
  std::stringstream tmp_ss5;
  tmp_ss5 << "./run"<< run <<"/gmm_models/UniversalGMM_Ng" << N_cent << "_dim" <<dim; 
  jo_gmm_diag universal_gmm_model;
  universal_gmm_model.load( tmp_ss5.str() );
  
  
  
  for (int pi = 0; pi <n_peo_test; ++pi ){ 
    
    std::stringstream tmp_ss4;
    tmp_ss4 << "./features/multi_test/feat_dim" << dim << "_" << peo_test(pi);  
    //tmp_ss4 << "./features/multi_test/feat_"<< peo_test(pi);  
    
    std::stringstream tmp_vec_lab;
    tmp_vec_lab << "./features/multi_test/lab_"<< peo_test(pi); 
    
    test_features_arma.load( tmp_ss4.str() );
    
    
    real_labels.load( tmp_vec_lab.str(), raw_ascii );
    
    
    int num_frames = test_features_arma.n_rows;
    
    cout << "num_frames: " << num_frames << endl;
    //getchar();
    
    mat log_prob;
    log_prob.zeros(num_frames, actions.n_rows); 
    
    cout << "Doing for person " << peo_test(pi) << endl;
    

    
    for ( int fr = 0; fr < num_frames -L +1 ; fr=fr+2)
    {
      //cout << fr << endl;
      mat mat_features;
      
      //cout << " j = " << fr << " to  j< " << fr+L << endl ;
      for (int j = fr; j< fr+L; j++)
      {
	//cout << j << " " ;
	mat_features	 = join_rows( mat_features, test_features_arma(j) );
	
      }
      
      //cout << "r&c "<<  mat_features.n_rows << " & " << mat_features.n_cols << endl;
      
      
      
      //cout << endl;
      mat likelihood_actions_lab = get_loglikelihoods(mat_features, run);
      likelihood_actions = likelihood_actions_lab.col(0);
      likelihood_label = likelihood_actions_lab.col(1);
      
      
      
      ///normalised log_likelihood:
//       universal_likelihood = universal_gmm_model.avg_log_p(mat_features);
//       vec vec_universal_likelihood;
//       vec_universal_likelihood.zeros(dim);
//       vec_universal_likelihood.fill(universal_likelihood);
//       vec nor_likelihood_actions  = likelihood_actions - vec_universal_likelihood;
//       log_prob.rows(fr, fr + L-1).each_row() += nor_likelihood_actions.t();
      //end nor_log_likelihood
      
      log_prob.rows(fr, fr + L-1).each_row() += likelihood_actions.t(); //comment this if you want to use the norm_log_likelihood

      
      
      //likelihood_actions_scenes.t().print();
      //cout << endl << "ini " << fr << " fin:" << fr + L-1 << endl;
      //
      
     
      
      
      //getchar();
      //cout << log_prob.rows(0, fr + L + 5) << endl;
      
      //mat X = log_prob.submat(0,0, fr + L + 5,6);
      //X.print();
      //X.submat( first_row, first_col, last_row, last_col )
      //getchar();
      
      
    }
    //log_prob.print();
    //cout << "Done here" << endl;
    //getchar();
    vec est_labels(num_frames);
    
    
    
    for ( int i = 0; i < num_frames; ++i)
    {
      //cout << " " << i;
      uword index;
      //double max_frame = 
      log_prob.row(i).max(index);
      //cout << "Log_prob_frame " << i << "= " << log_prob.row(i) << endl;
      
      //cout << "Max for " << log_prob.row(i) << " is " << max_frame << " pos: " << index<< endl;
      
      est_labels(i) = likelihood_label(index);
      //cout << "est_labels(i) " << est_labels(i) << endl;
      
      
      
      //getchar();
      
    }
    
    std::stringstream tmp_save_estlab;
    tmp_save_estlab<< "./run"<< run <<"/results/Ng_" << N_cent << "est_"<<peo_test(pi) <<"_dim" << dim << ".dat";  
    
    est_labels.save(tmp_save_estlab.str(), raw_ascii);
    
    
    std::stringstream tmp_save_real_lab;
    tmp_save_real_lab << "./run"<< run <<"/results/real_"<<peo_test(pi) <<".dat";  
    
    real_labels.save( tmp_save_real_lab.str(), raw_ascii );
    
    uvec comparing = find( est_labels == real_labels);
    double acc = comparing.n_elem;
    cout <<"performance for person "<<peo_test(pi)<<" is "<<setprecision(2)<<fixed<<100*acc/est_labels.n_elem << " %"<<endl;
    //Average performance not used as I'm running using 4 cores
    ave_acc = ave_acc + 100*acc/est_labels.n_elem;
    //getchar();
    
  }
  
  //cout << "Average performance is " << setprecision(2) << fixed << ave_acc/person.n_rows << " %"<<endl;
  
}


inline
mat
opt_feat::get_loglikelihoods(mat mat_features,  int run)
{
  
  //cout << "r&c "<<  mat_features.n_rows << " & " << mat_features.n_cols << endl;
  
  //cout << "get_loglikelihoods" << endl;
  mat likelihood_actions(actions.n_rows,2);
  vec likelihood_actions_lab(actions.n_rows); 
  vec likelihood_label(actions.n_rows);
  
  
  int pos = 0;
  
  
  for (uword act = 0; act < actions.n_rows; ++act)
  {
    jo_gmm_diag gmm_model;
    std::stringstream tmp_ss5;
    tmp_ss5 << "./run"<< run <<"/gmm_models/Ng" << N_cent << "_" << actions(act)<< "_dim" <<dim; 
    //tmp_ss5 << "./run"<< run <<"/gmm_models/Ng" << N_cent << "_" << actions(act); 
    //cout << tmp_ss5.str() << endl;
    gmm_model.load( tmp_ss5.str());
    
    likelihood_actions_lab (pos) = gmm_model.avg_log_p(mat_features);
    likelihood_label (pos) = act;
    pos++;
  }
  
  
  likelihood_actions.col(0) = likelihood_actions_lab;
  likelihood_actions.col(1) = likelihood_label;
  
  //likelihood_actions.print();
  //cout << "Done. Press a key" << endl;
  //getchar();
  return likelihood_actions;
  
}

